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1.  SUMMARY 


With  the  rapid  increase  in  seismic  data  collected  around  the  world,  there  is  an  acute 
need  to  integrate  the  available  data  to  construct  an  accurate  and  comprehensive  Earth 
model.  In  this  project,  we  develop  a  novel  framework  to  construct  Earth  models  by 
combining  the  multigrid/multilevel  method  and  full-wave  tomography:  the  multigrid 
method  is  designed  for  numerical  efficiency,  while  full-wave  tomography  is  well  suited 
for  iterative  update  of  the  3D  model  and  integration  of  various  seismic  data.  This  new 
framework  provides  an  efficient  way  to  progressively  improve  the  Earth  model  to  as  high 
resolution  as  observational  constraints  warrant.  We  apply  this  framework  to  build  a  3D 
velocity  model  for  the  eastern  hemisphere  that  covers  Eurasia,  Africa,  the  Middle  East, 
Australia,  and  the  Indian  Ocean  (latitude  55°S-55°N;  longitude  30°W-158°E). 

Using  a  frequency-time  data  normalization  procedure,  we  show  that  it  is  possible  to 
extract  up  to  300-s  period  Rayleigh  waves  from  ambient  noise  recorded  by  temporary 
broadband  seismic  stations  and  up  to  600-s  period  signals  from  permanent  stations  with 
very  broadband  sensors.  These  very  broadband  waves  extend  surface  wave  constraints  to 
all  depths  of  the  crust  and  upper  mantle.  We  collected  and  processed  up  to  23  years 
(1990-2012)  of  continuous  seismic  data  from  all  available  broadband  seismic  stations  in 
the  eastern  hemisphere.  Empirical  Green’s  functions  (EGEs)  extracted  from  ambient 
seismic  noise  are  used  in  full-wave  ambient  noise  tomography  developed  under  this 
project.  Three  levels  of  finite-difference  computation  and  inversion  have  been  carried 
out  at  the  hemispherical,  continental,  and  regional  scales  to  progressively  construct  a 
hierarchical,  multi-resolution  P  and  S  velocity  model.  Solutions  from  the  coarse  grids  are 
used  as  initial  models  for  the  finer  grids  and  high-resolution  models  from  the  finer  grids 
are  passed  to  the  coarse  grids  to  be  integrated  in  the  next  model  interation.  The  results 
show  much-improved  resolution  compared  to  that  of  conventional  ambient  noise 
tomography  and  earthquake  tomography.  The  new  model  provides  constraints  on  several 
geological  processes  in  the  region,  though  detailed  geological  interpretation  is  beyond 
this  report.  By  comparing  observed  and  synthetic  waveforms  calculated  from  selected 
velocity  models,  we  establish  a  procedure  to  directly  and  systematically  validate  the 
models.  Our  results  suggest  that  all  future  models  should  undergo  rigorous  model 
validation  using  independent  seismic  observations. 

We  extend  the  3D  Green’s  function  based  source  moment  inversion  from  local  to 
regional  scales  and  obtain  moment-tensor  determinations  of  selected  earthquakes  in  the 
southeast  Tibetan  plateau.  The  solutions  are  in  general  agreement  with  the  global  centroid 
moment  tensor  solutions,  though  the  two  source  locations  may  differ  by  10s  of  km. 
Inspection  of  waveform  fit  suggests  that  regional  stations,  the  3D  reference  model  and 
short-period  body  waves  provide  tighter  constraints  on  the  locations  of  the  earthquakes. 
The  lowest  magnitude  of  the  earthquakes  in  this  study  is  Mw  4.35.  As  the  crustal  velocity 
model  improves,  short-period  (<15  s)  surface  waves  become  useful  and  the  SGT-based 
method  may  be  applied  to  smaller  earthquakes  at  local  and  regional  scales. 
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2.  INTRODUCTION 


There  is  a  general  eonsensus  that  the  use  of  3D  Earth  models  improves  seismic 
characterization  (Antoun  et  ah,  2008;  Covellone  and  Savage,  2012).  Advance  in  high- 
performance  computation  has  made  it  feasible  to  simulate  propagation  of  increasingly 
broadband  seismic  wave  in  3D  models.  Nevertheless  it  remains  impractical  to  routinely 
compute  forward  wave  propagation  from  seismic  sources  in  3D  regional  models  in  real 
time.  One  way  to  overcome  this  challenge  is  to  take  advantage  of  an  important  property 
in  seismic  wave  propagation,  the  source-receiver  reciprocity.  The  use  of  the  reciprocity 
property  and  Green’s  functions  can  drastically  reduce  computational  costs  when  the 
sources  outnumber  the  stations  (e.g.,  Eisner  and  Clayton,  2001;  Graves  and  Wald,  2001). 
More  importantly,  with  3D  Green’s  functions  pre-calculated  and  saved  in  a  database,  they 
can  be  extracted  quickly  from  the  database  after  a  seismic  event  and  used  to  determine 
seismic  moments  (Zhao  et  ah,  2006;  Shen  et  ah,  2014). 

The  usefulness  of  3D  Green’s  functions  depends  on  the  accuracy  of  the  corresponding 
3D  Earth  models.  By  more  accurately  acounting  for  the  physics  of  wave  propagation  in 
the  3D  Earth,  the  full-wave  methods  have  been  proven  to  be  effective  in  improving 
model  resolution  (e.g.,  Chen  et  ah,  2007b;  Tape  et  ah,  2009;  Gao  and  Shen,  2014).  Two 
methods  have  been  developed  to  date  to  carry  out  full-wave  tomography.  One  is  the 
adjoint- wavefield  (AW)  method,  which  back-propagates  data  from  receivers  to  image 
structure  (e.g.,  Tromp  et  ah,  2005;  Eiu  and  Tromp,  2006).  The  other  is  the  strain  Green’s 
tensor  (SGT)  based,  scattering-integral  (SI)  method,  which  calculates  and  stores  the 
sensitivity  kernels  for  each  data  functional  (e.g.,  Zhao  et  al,  2005;  Chen  et  ah,  2007a,b; 
Zhang  et  ah,  2007).  Both  methods  are  based  on  simulation  of  full-wave  propagation  in 
3D  models.  The  main  differences  are  in  computation.  In  general,  the  AW  is  CPU 
intensive,  while  the  SI  approach  requires  a  large  disk  space  and  fast  input/output.  One 
advantage  of  the  SI  approach  is  that  the  SGT  database  can  be  used  for  source  moment 
inversion  in  an  efficient,  global  optimization  scheme  (Shen  et  ah,  2014).  Eor  this  reason, 
we  focus  on  the  development  and  application  of  the  SGT-based  approach  in  this  project. 

Eigure  1  outlines  the  workflow  of  the  unified  tomography  and  source  moment 
inversion  based  on  3D  SGT  databases.  Station  SGTs  are  constructed  from  a  3D  reference 
model  by  finite-difference  simulation  of  the  responses  to  orthogonal  unit  impulsive  point 
forces  acting  at  the  stations.  By  extracting  the  SGTs  in  a  small  volume  surrounding  the 
source  reference  location,  we  can  invert  for  source  moment  tensors  and  location  in  a 
global  optimization  scheme  (Shen  et  ah,  2014).  The  same  SGTs  can  also  be  used  to 
calculate  the  sensitivities  of  waveforms  and  their  derivatives  (phase  delays  and  amplitude 
anomalies)  to  source  parameters.  Eorward  wave  propagation  simulation  from  the  source 
in  the  3D  reference  model  provides  synthetic  waveforms  at  stations  and  the  wavefield 
within  the  volume  of  the  model.  Travel  time,  amplitude,  or  waveform  anomalies  are 
measured  from  the  observed  and  synthetic  waveforms  at  stations.  The  forward  wavefield 
(from  the  source)  and  the  station  SGTs  are  used  to  calculate  the  finite-frequency 
sensitivities  to  perturbations  in  Vp,  Vs  (or  bulk  and  shear  moduli),  density,  and 
attenuation.  Together  with  the  sensitivities  to  source  parameters,  the  measurements  and 
structural  sensitivity  kernels  are  used  to  invert  for  the  earth  structure  (e.g.,  Zhang  et  ah, 
2009).  The  tomographic  inversion  results  and  additional  constraints  (e.g.,  receiver 
function  solutions,  well  logs,  etc)  are  added  to  the  3D  reference  model.  This  process  can 
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be  repeated  to  progressively  improve  the  solution.  For  full-wave  tomography  using 
ambient  noise  data,  the  souree  moment  inversion  is  omitted  and  the  synthetie  waveforms 
at  stations  and  sensitivity  kernels  are  all  ealeulated  from  the  SGTs  (Shen  and  Zhang, 
2010;  Gao  and  Shen,  2014). 


3D  structural  models 


source  kernels 


Figure  1.  Flow  chart  of  the  unified  tomography  and  moment  tensor  inversions  based  on  3D  finite- 
difference  SGT  databases. 
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One  of  the  main  challenges  in  full-wave  tomography  is  the  intense  computation 
associated  with  simulating  wave  propagation  in  3D  models.  To  reduce  computational 
cost,  which  generally  scales  with  the  number  of  discrete  nodes  in  numerical  methods,  we 
adapt  a  multigrid/multilevel  method  (Briggs, 

1987)  to  solve  the  wave  propagation  and 
inversion  problem  (Figure  2).  We  use 
several  levels  of  finite-difference 
computation  with  increasingly  finer  grids  at 
the  hemispherical,  continental,  regional,  and 
local  scales  to  progressively  construct  a 
hierarchical,  multi-resolution  P  and  S 
velocity  model  (Figure  3).  Data  analysis 
correspondingly  evolves  from  long-period 
waves  for  teleseismic  observations  to 
broadband  waves  for  regional  and  local 
observations.  The  solution  of  a  lower 
(coarser)  grid  provides  the  starting  model  for 
the  next  higher  (finer)  level.  The  finer-scale 
solution  is  then  passed  to  the  coarse  model 
for  the  next  model  iteration  in  an  iterative. 


number  of  iteration 

Figure  2.  Schematic  illustration  of  the  multigrid 
full-wave  tomographic  framework.  Inversion 
starts  from  a  coarse  grid.  When  the  solution 
converges  after  several  iterations  (usually  3  to 
4  in  our  experience),  computation  is  moved  to  a 
finer  grid  in  a  smaller  area.  The  solution  is 
passed  from  coarse  to  fine  grids,  and  vice 
versa. 
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multi-grid  scheme  (Briggs,  1987).  Within  this  framework,  the  levels  of  grids  and 
eorresponding  multigrid  resolution  ean  vary  from  the  global  to  any  survey  seales,  limited 
only  by  the  availability  of  seismie  data  and  eomputational  resourees. 

The  iterative  update  of  the  3D  Earth  model  in  full-wave  tomography  (Figure  1)  is 
well  suited  for  the  multigrid  method.  Together,  they  provide  an  effieient  way  to 
eumulatively  and  progressively  improve  the  Earth  model.  Although  eaeh  of  the  lower- 
level  models  in  this  study  eneompasses  one  and  only  one  next-level  model  (Figure  3), 
several  branehes  of  the  multigrid  model  ean  be  developed  in  parallel  (e.g.,  in  Afriea,  the 
Middle  East,  and  Australia). 


Longitude 

Figure  3.  A  multigrid  approach  to  develop  a  joint  P-  and  S-velocity  model  in  the  eastern 
hemisphere.  The  three  levels  of  grids  represent  the  hemispherical-,  continental-,  and  regional- 
scale  models.  Red  triangles  are  selected  permanent  and  temporary  broadband  seismic  stations 
used  in  inversion.  Grey  lines  are  the  coast  and  the  2000-m  and  4000-m  topographic  contours. 

The  products  of  this  project  include:  (1)  a  new  non-staggered-grid  finite-difference 
code  that  simulates  wave  propagation  in  a  heterogeneous  and  anisotropic  Earth  in 
spherical  coordinates  (Zhang  et  ah,  2012),  (2)  a  frequency-time  normalization  method 
that  significantly  improves  the  quality  of  the  empirical  Green’s  functions  (EGFs) 
extracted  from  ambient  seismic  noise  (Shen  et  ah,  2012)  and  provides  very  broadband 
signals  that  enable  imaging  of  the  entire  crust  and  upper  mantle,  (3)  a  full-wave  ambient 
noise  tomographic  model  of  the  eastern  hemisphere,  (4)  a  procedure  to  validate  velocity 
models  (Gao  and  Shen,  2012;  Gao  and  Shen,  in  preparation),  and  (5)  SGT-based  moment 
tensor  inversion  (Shen  et  ah,  2014).  In  the  following,  we  focus  on  the  new  ambient  noise 
data  processing  method  (section  3),  full-wave  ambient  noise  tomography  in  the  eastern 
hemisphere  (section  4),  model  validation  (section  5),  and  SGT-based  moment  inversion 
(section  6).  Concluding  remarks  and  recommendations  are  presented  in  Sections  7  and  8, 
respectively. 
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3.  EXTRACTION  OF  VERY  BROADBAND  EMPIRICAL  GREEN’S 
FUNCTIONS 

Advance  in  seismic  imaging  comes  from  two  main  sourees:  new  data  sets  and  new 
inverse  methods.  The  extraetion  of  EGFs  from  ambient  seismic  records  is  a  good 
example  of  new  data  sets:  EGFs  provide  powerful  eonstraints  on  the  Earth’ s  strueture, 
partieularly  in  aseismie  regions,  beeause  they  do  not  depend  on  earthquake  sources,  have 
preeisely  known  virtual  souree  loeations,  and  yield  useful,  shorter  period  surfaee  waves 
than  from  earthquakes.  Previous  ambient  noise  tomographie  studies  use  short-  to 
intermediate-period  (up  to  '-50  s)  surfaee  waves  that  provide  constraints  mainly  in  the 
crust  and  uppermost  mantle  (e.g.,  Yao  et  ah,  2006;  Ein  et  ah,  2008;  Yang  et  ah,  2008).  In 
this  projeet,  we  develop  a  new  data  proeessing  method  to  extraet  very  broadband 
Rayleigh  waves  from  ambient  seismie  reeords,  extending  EGF  eonstraints  to  the  entire 
erust  and  upper  mantle  (Shen  et  ah,  2012;  Shen  et  ah,  2013). 

Theoretieal  and  experimental  studies  have  shown  that  retrieval  of  Green’s  funetions  is 
best  aehieved  in  a  diffuse  wavefield,  in  which  amplitudes  of  normal  modes  are 
uneorrelated  and  have  equal  power  (Weaver  &  Lobkis  2001)  or  energies  radiated  by 
random  sourees  are  at  the  same  level  (e.g.,  Snieder,  2004;  Curtis  &  Holliday  2010). 
However,  original  seismie  reeords  are  highly  non-stationary  and  eontain  local  noise, 
instrument  irregularities,  and  earthquake  arrivals.  To  eompensate  for  the  deviation  from 
equal  partitioning  of  noise  sourees,  Stehly  et  al.  (2008)  showed  that  the  use  of  the  eoda  of 
eross-eorrelation  funetions  helps  remediate  the  problem  of  anisotropie  distribution  of 
noise  sourees.  For  applications  with  a  regular  and  well-sampled  array  of  stations, 
multidimensional  deeonvolution  (Wapenaar  et  al,  2011)  and  directional  balancing 
(Curtis  &  Halliday,  2010)  provide  alternatives  to  the  eross-eorrelation  method.  But  for 
many  passive  seismie  interferometry  problems,  eross  eorrelation  of  reeords  at  station 
pairs  remains  a  eommonly  used  approach. 

Bensen  et  al.  (2007)  provided  an  excellent  review  of  the  cross-eorrelation-based 
methods.  The  two  most  representative  practiees  of  data  normalization  are  one-bit 
normalization,  in  which  the  time  series  is  normalized  by  its  point-by-point  absolute  value 
(e.g.,  Shapiro  et  al,  2005),  and  running-absolute-mean  normalization,  in  whieh  the  reeord 
is  normalized  by  a  running  average  of  its  windowed  absolute  value  (Bensen  et  al,  2007). 
Although  long-period  surface  waves  have  been  observed  with  these  methods  (e.g., 

Shapiro  &  Campillo  2004;  Bensen  et  al,  2007;  Nishida  et  al.,  2009),  ambient  noise 
studies  using  data  from  temporary  broadband  experiments  typieally  have  useful  Rayleigh 
waves  at  periods  up  to  50  s  (e.g.,  Yao  et  al.,  2006;  Yang  et  al.,  2008).  Using  a  frequeney- 
time  normalization  proeedure  developed  under  this  project,  we  show  that  mueh  broader 
frequeney  Rayleigh  waves  ean  be  extraeted  (Shen  et  al,  2012). 

We  have  eolleeted  and  processed  continuous  data  recorded  by  broadband  seismie 
stations  in  the  eastern  hemisphere  between  1990  to  2012  (Shen  et  ah,  2013).  Data 
proeessing  follows  a  nested  approaeh  that  eorresponds  to  the  multi-levels  of  the  study 
areas  (Figure  3,  Table  1).  For  the  hemispherieal-seale  model,  for  example,  we  seleeted 
'-'290  stations  having  at  least  300  days  of  eontinuous  reeords  for  further  analysis.  Figure  4 
shows  an  example  of  EGFs  from  eross  eorrelation  of  vertieal-vertieal  eomponents  of 
permanent  seismie  stations.  Our  results  eonfirm  that  elear  Rayleigh  wave  signals  ean  be 
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extracted  from  seismic  “hum”,  or  earth’s  background  free  oscillation  (Nishida  et  ah, 
2009).  Furthermore  we  extend  the  useful  period  from  400  s  (Nishida  et  ah,  2009)  up  to 
600  s!  The  remarkably  good  symmetry  of  very-long  period  (>100  s)  EGFs  for  the 
positive  and  negative  time  lags  of  cross  correlation  indicates  a  relative  homogeneity  of 
noise  sources  at  low  frequencies.  For  PASSCAL  broadband  stations  with  1-2  years’ 
records,  we  can  extract  useful  Rayleigh  waves  up  to  300  s  period  (Shen  et  ah,  2012). 


Figure  4.  Example  of  EGFs  derived  from  ambient  noise  (Earth  “hum”)  with  various  wave 
periods  (upper  left:  50-600  s;  upper  right:  100-600  s;  lower  left:  200-600  s;  and  lower  right:  300- 
600  s).  The  lines  in  the  accompanying  maps  show  the  paths  with  high  signal-to-noise  arrivals 
from  the  virtual  source  (station  ElYB)  to  other  receivers  in  the  eastern  hemisphere.  The  red  dots 
mark  a  wave  propagating  at  a  velocity  of  3.9  km/s. 


Table  1.  Model  Information 


Grid/Model  Level 

number  of 

seismic 

stations 

EGF  period 
range 

number  of 
phase 

measurements 

Lateral  grid 
spacing  (deg.) 

Inversion  cell 
lateral  spacing 
fdeg.) 

Eastern  hemisphere 

290 

50-600  s 

23,000 

0.2 

2.0 

Western  China 

415 

25-300  s 

11,000 

0.1 

1.0 

Tibetan  plateau 

400 

12.5-100  s 

14,000 

0.025 

0.5 
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4.  FULL-WAVE  TOMOGRAPHY 


4,1  Reference  Model  and  Wave  Propagation  Simulation 

The  initial  3D  reference  model  for  the  lowest-level  grid  is  composed  of  CRUST  2.0 
(Bassin  et  ah,  2000)  for  the  crust  and  AK135  (Kennett  et  ah,  1995)  for  the  mantle.  The 
choice  of  an  average  Earth  model  for  the  mantle  instead  of  one  of  the  current  3D  models 
(e.g.,  Ritsema  and  van  Heijst,  2000;  Ritzwoller  et  ah,  2002;  Antolik  et  ah,  2003)  is  a 
tradeoff  between  finding  a  good  starting  model 
and  the  need  to  explore  a  wide  enough  model 
space  and  verify  the  final  solution. 

We  use  a  non-staggered-grid  finite-difference 
code,  a  code  developed  under  this  project,  to 
simulate  wave  propagation  in  spherical 
coordinates  (Figure  5,  Zhang  et  ah,  2012).  The 
DRP/opt  MacCormack  scheme,  which  alternately 
uses  forward  and  backward  finite-difference 
operators,  is  implemented  in  the  calculation.  The 
solution  is  optimized  for  both  dispersion  error  and 
dissipation  error,  and  has  4th-order  dispersion 
accuracy.  In  order  to  increase  the  efficiency  of 
the  finite-difference  algorithm,  we  use  a  grid  with 
non-uniform  grid  spacing  to  discretize  the 
computational  domain.  The  grid  varies 
continuously  with  smaller  spacing  in  low  velocity 
layers  and  with  larger  spacing  otherwise.  The 
computational  domain  is  surrounded  by  complex- 
frequency  shifted  Perfect  Matched  Layers  (PMLs) 
implemented  through  auxiliary  differential 
equations  (Zhang  and  Shen,  2010). 


Figure  5.  Geometry  of  typical  non- 
uniform,  non-staggered  finite -difference 
grid  in  a  spherical  coordinate  system. 


For  the  lowest-level  grid,  the 
computational  domain  covers  the  eastern 
hemisphere  (latitude  55°S-55°N; 
longitude  30°W-158°E)  with  a  lateral 
grid  spacing  of  0.2°,  and  a  non-uniform 
grid  from  the  surface  to  --2100  km  depth 
(3  km  near  surface,  gradually  increasing 
to  ~40  km  near  the  bottom  of  the 
model).  Tests  indicate  that  the  numerical 
dispersion  errors  are  less  than  0.1%  for 
Rayleigh  waves  with  periods  greater 
than  50  s.  Figure  6  shows  a  snapshot  of 
simulated  seismic  response  to  a  vertical 
force  at  station  BJI  (Beijing).  The  two 
higher-level  grids  (levels  2  and  3)  are 
spaced  at  0.1°  and  0.025°,  respectively. 


Figure  6.  A  snapshot  of  the  seismic  response  at  the 
surface  (vertical  ground-motion  velocity)  to  a 
vertical  force  at  station  BIT  (Beijing). 
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with  correspondingly  shorter  usable  periods  of  synthetic  waveforms  (Table  1). 

4,2  Phase  Delay  Measurements 

We  measure  Rayleigh  wave  phase  delays  by  cross-correlating  observed  and  synthetic 
Green’s  functions.  For  the  hemispherical-scale  model,  the  waveforms  are  band-pass 
fdtered  in  six  period  ranges  (50-100,  75-150,  100-200,  150-300,  200-400,  and  300-600  s). 
A  total  of -23,000  phase  delay  measurements  meet  our  data  seleetion  criteria  (waveform 
correlation  coefficient  greater  than  0.7  and  signal-to-noise  ratio  greater  than  5).  Most  of 
the  propagation  paths  with  useful  waves  are  in  the  northern  hemisphere  (Figure  7). 
Nevertheless  the  expansion  of  the  initial  study  area  to  include  stations  in  the  southern 
hemisphere  provides  the  north-south  propagation  paths  that  are  essential  to  obtain  good 
resolution  in  the  north.  Table  1  lists  the  number  of  stations,  frequency  range,  and  number 
of  phase  measurements  for  other  grid  levels. 


0  50  100  150 

Long. 


0  50  100  150 

Long. 


Figure  7.  Great  circle  paths  of  useful  Rayleigh  wave  delays  from  EGFs  at  200-400  s  (left), 
and  75-150  s  (right)  periods.  Triangles  mark  the  locations  of  seismic  stations.  Eurasia, 
northern  Afiica  and  Australia  are  well  sampled  at  all  periods  used. 

4,3  The  3D  Sensitivity  Kernels 

The  cross-correlation  of  vertical-vertical  components  at  pairs  of  seismic  stations 
yields  predominantly  Rayleigh  waves  in  the  time  window  and  frequency  of  our  interest. 
Although  commonly  neglected  in  surface  wave  tomography,  P-wave  speed  affeets 
Rayleigh  waves  (Figure  8).  The  sensitivity  to  P-wave  speed  is  near  the  surfaee  (the  depth 
of  sensitivity  increases  with  period)  and  at  a  level  comparable  to  the  Vs  kernel  near  the 
surface.  The  sensitivity  kernels  are  calculated  with  the  SGT-based,  scattering-integral 
method  (Zhao  et  ah,  2005;  Chen  et  ah,  2007a;  Zhang  et  ah,  2007). 


Va.  hcriontBl  al  253  M  km  4eptti  V*.  honzcrW  Ml  at  2S3.S9  lim 


Figure  8.  Finite-frequency  sensitivity  of  long-period  EGFs  between  BIT  (Beijing)  and  LSA  (Lhasa)  to  Vp  (left) 
and  Vs  (right)  perturbations.  The  Rayleigh  wave  period  is  150-300  s.  The  horizontal  slice  is  at  ~250  km  depth. 
The  sensitivity  to  Vp  is  mostly  in  the  upper  100  km.  The  unit  of  the  sensitivity  kernel  is  s/m^. 
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4,4  Inversion  and  the  Velocity  Model 

We  represent  the  Rayleigh  wave  delays  as  a  joint  Vp  and  Vs  inverse  problem 
linearized  relative  to  an  iteratively  updated  3D  model, 

(5t  =  /[R„(mo,x)Am„+/<^i(mo,x)Am/j]a'\/  (1) 

where  5t  is  the  travel  time  delay  relative  to  the  synthetics  for  the  3D  model  mo,  Ama  and 
Amp  perturbations  to  P-wave  and  S-wave  speeds,  respectively,  and  the  sensitivity 
kernels  to  P-wave  and  S-wave  speeds,  respectively.  The  integration  is  for  the  3D  volume 
of  the  model.  The  sensitivity  of  Rayleigh  waves  to  density  is  not  explicitly  expressed  in 
the  above  equation,  but  density  is  updated  after  each  model  iteration  using  a  density-Vp 
relation  (Christensen  and  Mooney,  1995). 

The  above  equation  is  discretized  into  inversion  cells.  For  the  hemispherical-scale 
inversion,  we  use  2.0°x2.0°  blocks 
as  the  inversion  cells.  The 
thickness  of  the  blocks  changes 
from  ~14  km  near  the  surface  to 
~80  km  at -1600  km  depth.  The 
inversion  block  size  is  reduced  for 
the  higher-level  grids  to  capture 
sharp  velocity  changes  (Table  1). 

The  size  of  the  inversion  cells  for 
level  3  (the  Tibetan  plateau)  is 
0.5°x0.5°.  The  thicknesses  of  the 
blocks  are  proportional  to  their 
lateral  dimensions. 

The  inverse  problem  is  solved 
with  LSQR  (Paige  and  Saunders, 

1982)  with  damping  and 
smoothness  constraints.  The 
damping  and  smoothness 
constraints  are  gradually  reduced 
until  the  data  misfit  is  consistent 
with  the  estimated  data  uncertainty 
(normalized  x  ~1),  as  in  Montelli 
et  al.  (2004).  The  uncertainty  of 
the  phase  delays  is  calculated  from 
the  monthly  EOF  stacks. 

Figure  9  is  a  comparison  of  the 
Vsv  model  from  this  study  and  that 
of  Shapiro  and  Ritzwoller  (2002), 
which  was  derived  from  group  and 
phase  velocities  of  earthquake  surface  wave  arrivals.  Despite  the  completely  different 
datasets  and  methods  used  in  the  two  studies,  there  are  remarkable  similarities  between 
the  two  models.  For  example,  the  African  cratons  and  the  west  Australian  craton  exhibit 


Vs  (This  study,  depth  =  198  km) 
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Figure  9.  Comparison  of  the  shear-wave  velocity  models  at 
198  km  depths  from  full-wave  ambient  noise  tomography 
(upper  panel)  and  CUB  2.0  (lower  panel,  Shapiro  and 
Ritzwoller,  2002).  Color  indicates  the  perturbation  in 
percentage  relative  to  the  average  at  the  depth.  Grey  and 
green  lines  mark  the  coast  and  plate  boundaries. 
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high  velocity  anomalies  in  both  models.  The  Afar  and  Eastern  Africa  Rift  are  prominent 
low  velocity  anomalies  in  both  models.  There  are  also  notable  differences  between  the 
two.  The  full-wave  ambient  noise  model  has  resolution  in  the  mantle  transition  zone 
(Figure  10),  while  resolution  in  CUB2.0  is  limited  to  the  upper  ~300  km  of  the  mantle 
due  to  the  limited  period  range  of  the  earthquake  arrivals  (<  150  s  period).  The  full-wave 
ambient  noise  model  also  reveals  far  more  detailed  structure  than  in  CUB2.0,  particularly 
in  west  China  and  the  Tibetan  plateau  (Figures  9  and  11),  where  full-wave  ambient  noise 
tomography  with  fine  grids  and  dense  PASSCAF  stations  yields  much  higher  resolution. 
Another  notable  major  difference  between  the  two  models  is  the  broad  and  strong  high 
velocity  anomaly  beneath  the  Tibetan  platan  in  CUB2.0.  In  our  new  model,  the  high 
velocity  anomalies  are  fragmented  beneath  the  plateau,  consistent  with  body- wave 
tomography  results  (Fiang  et  ah,  2012),  and  more  prominent  in  the  areas  surrounding  the 
plateau,  as  expected  from  the  presence  of  the  Sichuan  Basin  (the  Yangtze  craton)  and  the 
Ordos  Basin  (the  North  China  craton)  east  of  the  plateau,  the  Tarim  Basin  (the  Tarim 
Archean  craton)  north  of  the  plateau,  and  the  Indian  plate  (craton)  in  the  south. 


Vs  (full  wave,  d  =  575  km) 


Figure  10.  Shear-wave  velocity  (Vsv)  at  575  km  depth  from  full-wave  ambient  noise  tomography. 
Color  indicates  the  perturbation  in  percentage  relative  to  the  average  at  the  depth. 
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Figure  1 1 .  (a)  Broadband  seismic  stations  (triangles)  used  in  full-wave  ambient  noise  tomography 
of  the  Tibetan  plateau.  The  3000-m  contour  outlines  the  plateau,  (b)  The  shear-wave  velocity  at 
depth  of  ~46  km.  The  color  scale  (red  to  blue)  is  between  3.59-4.59  km/s.  Dashed  white  lines  are 
the  sutures  in  Tibet,  while  dashed  grey  lines  are  faults.  The  velocity  structures  closely  follow  the 
geological  features.  As  an  indication  of  resolution,  the  ductile  extmsion  of  high-grade 
metamorphic  rocks  on  the  southern  slope  of  the  Himalaya  is  imaged  as  a  band  of  high  velocity 
anomaly. 
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Synthetic  inversion  of  various  input  models  is  a  common  practice  in  tomographic 
resolution  analysis.  This  approach  has  limitations  and  can  be  misleading  in  the  sense  that 
synthetic  inversion  explores  only  a  limited  model  subspace  (Leveque  et  ah,  1993). 
Nevertheless,  it  is  useful  if  interpretations  are  restricted  to  the  model  subspaces  explored. 
With  this  caveat  in  mind,  we  run  3D  checkerboard  resolution  tests  with  a  maximum  of 


±5%  velocity  perturbations  for  both  Vp 
and  Vs  at  the  various  grid  levels.  The 
synthetic  phase  delay  times  are  calculated 
with  the  sensitivity  kernels  (Eq.  1).  The 
uncertainties  of  individual  observed  phase 
measurements  are  estimated  from  the 
variations  of  monthly-stacked  EGFs  and 
added  to  the  synthetic  phase  delays.  Figure 
12  shows  that  the  results  of  a  checkerboard 
test  at  the  hemisphere  level.  The  6°x6° 
checkerboard  velocity  perturbations  are 
well  reconstructed  in  Eurasia,  northern 
Africa,  Australia,  and  the  northern  Indian 
Ocean.  For  4°x4°  checkerboard  velocity 
perturbations  (not  shown),  western  China, 
the  Middle  East,  the  Mediterrainean,  North 
Africa,  Australia,  and  the  northern  Indian 
Ocean  remain  well  resolved,  though 
smearing  becomes  obvious  at  the  margins. 
For  the  reginal  scale  (level  3), 
checkerboard  tests  show  that  l°xl° 
velocity  pertubations  are  well  resolved  in 
eastern  and  southern  Tibet  (not  shown). 


Vs  at  depth  •  163.2  km 
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Figure  12.  Resolution  at  the  hemisphere  level. 
The  input  (lower  panel)  and  recovered  (upper 
panel)  checkerboard  velocity  perturbations  at  1 63 
km  depth  show  good  resolution  in  Eurasia, 
northern  Africa,  Australia,  and  the  northern 
Indian  Ocean. 


5.  MODEL  VERIFICATION 
AND  VALIDATION 

There  has  been  an  increasing  number  of  published  seismic  velocity  models  each  year. 
With  models  resolved  from  different  methods  and  datasets  comes  the  need  to  understand 
the  accuracy  and  resolution  of  these  models.  Becker  (2012)  quantitatively  documented 
the  similarities  and  differences  among  several  velocity  perturbation  models  resolved  from 
body-wave  tomography  in  the  Pacific  Northwest.  Although  the  large-scale  structures  are 
found  to  be  consistent  among  these  models,  the  magnitudes  of  the  velocity  anomalies 
vary  within  a  wide  range,  posing  a  problem  in  using  the  velocity  models  to  infer  the 
structure  and  dynamics  of  the  mantle  in  the  region. 

Although  the  comparison  of  various  velocity  models  is  helpful,  it  does  not  provide  an 
assessment  of  the  validity  of  the  models  in  predicting  seismic  observations.  A  more 
straightforward  and  powerful  way  to  validate  models  is  to  directly  compare  the  observed 
and  synthetic  waveforms  (e.g..  Song  and  Helmberger,  2007;  Qin  et  ah,  2009;  Bozdag  and 
Trampert,  2010;  Gao  and  Shen,  2012).  Using  simulation  of  full  wave  propagation,  Gao 
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and  Shen  (2012)  systematically  tested  four  surface-wave  tomographic  models  in  the 
Cascades.  The  results  showed  that  none  of  the  Cascade  models  is  adequately  accurate  in 
the  sense  that  the  misfits  between  the  observed  and  predicted  seismic  waves  are  much 
larger  than  the  uncertainty  of  the  observed  waveform. 

Figure  13  shows  a  comparison  of  the  phase  delays  as  a  function  of  station  separation 
distance  at  several  frequency  bands  for  CUB2.0  (Shapiro  and  Ritzwoller,  2002)  and  our 
new  model  of  the  Tibetan  plateau.  At  long  periods  (50-100  s),  phase  delays  for  CUB2.0 
are  generally  positive  and  have  a  positive  trend  with  station  separation,  indicating  that 
the  model  is  too  fast  for  the  real  structure  beneath  the  Tibetan  plateau.  Our  new  model 
predicts  smaller  phase  delays  that  are  distributed  closely  around  zero.  In  order  to  further 
understand  full-wave  tomography,  Gao  and  Shen  (manuscript  in  preparation)  carried  out 
a  systematic  validation  of  several  selected  models  in  north  America,  where  the  US  Array 
and  unrestricted  access  to  common  data  sets  make  the  comparison  relatively 
straightforward.  The  results  show  that  the  full-wave  tomographic  model  provides  better 
prediction  for  both  earthquake  and  EGF  arrivals. 


50-100  s  35-70  s 


0  500  1000  1500  2000  0  500  1000  1500  2000 

Distance,  km  Distance,  km 

Figure  13.  Phase  delays  between  the  EGFs  from  the  stations  in  the  Tibetan  plateau  and  the 
synthetie  Green’s  functions  for  CUB2.0  (black  dots)  and  our  new  model  (red  dots)  in 
several  period  ranges.  The  error  bars  represent  1  standard  deviation  of  measurement  errors. 
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6.  SGT-BASED  MOMENT  TENSOR  INVERSION 

6.1  Overview  of  Moment  Tensor  Solutions  Based  on  3D  Green’s  Functions 

Advances  in  computation  have  made  the  routine  use  of  3D  Green  functions  in 
earthquake  source  inversion  an  increasingly  attractive  option.  Synthetic  waveforms  from 
good  3D  reference  models  account  for  complex  waveform  phenomena  and  reduce  phase 
differences.  Liu  et  al.  (2004)  determined  the  source  mechanisms  of  small  to  moderate 
earthquakes  in  southern  California  using  synthetic  waveforms  calculated  with  the 
spectral-element  method.  Because  the  derivatives  of  the  source  parameters  are 
determined  numerically  by  differentiating  synthetics  with  respect  to  the  source 
parameters,  up  to  10  forward  simulations  are  needed  for  each  earthquake  (6  moment 
tensor  components,  latitude,  longitude,  depth,  and  the  reference  location).  When  the 
reference  location  is  far  from  the  true  location,  the  derivatives  of  the  location  parameters 
may  not  adequately  reflect  the  non-linear  variation  of  waveforms  as  a  function  of 
location. 

Taking  advantage  of  source-receiver  reciprocity,  Zhao  et  al.  (2006)  introduced  the  use 
of  strain  Green  tensors  (SGTs)  in  source  inversion  based  on  3D  reference  models.  In  this 
approach,  the  SGTs  are  calculated  with  a  finite-difference  code  and  saved  in  a  database. 
The  number  of  forward  waveform  simulations  is  three  times  the  number  of  stations  used 
in  source  inversion.  When  the  sources  outnumber  the  stations,  the  Green’s  tensor  or  SGT- 
database  approach  is  more  computationally  efficient  (e.g..  Graves  &  Wald,  2001).  More 
importantly  since  the  SGTs  are  pre-calculated,  the  SGT-database  approach  is  better  suited 
for  (near-)  real-time  source  inversion.  Zhao  et  al.  (2006)  adopted  a  perturbation  approach 
in  source  inversion  and  chose  to  minimize  frequency-dependent  arrival  time  and 
amplitude.  In  this  project,  we  extend  the  SGT-database  approach  from  local  scales  with 
waveforms  calculated  in  Cartesian  coordinates  (Zhao  et  ah,  2006)  to  regional  and 
teleseismic  scales  with  waveforms  calculated  in  spherical  coordinates  (Zhang  et  ah, 
2012),  so  we  are  able  to  use  any  combination  of  local  and  teleseismic  waveforms  in  a 
source  inversion  (Shen  et  ah,  2014).  Instead  of  frequency-dependent  arrival  time  and 
amplitude  anomalies,  we  fit  multi-frequency  filtered  and  time-shifted  waveforms  directly. 
For  a  given  source  location,  the  waveforms  are  linear  functions  of  the  moment  tensor 
elements,  which  can  be  solved  as  a  linear  inverse  problem.  To  account  for  perturbations 
in  earthquake  source  location,  we  carry  out  a  grid  search  in  the  vicinity  of  the  reference 
location.  This  is  necessary  for  earthquakes  with  location  errors  of  5  km  and  more,  as  is 
the  case  for  most  earthquake  catalogs.  When  the  true  source  location  is  far  from  the 
reference  location,  the  derivative  of  waveform  misfit  about  the  reference  location  may 
not  be  a  good  approximation  to  obtain  the  true  location.  This  in  turn  affects  the 
waveforms  and  the  moment  tensor  solution.  The  grid  search  is  computationally  efficient 
as  it  involves  extracting  pre-calculated  SGTs  in  a  small  volume  containing  the  reference 
source  location  from  the  database.  We  apply  the  new  method  to  southeast  Tibet  in 
conjunction  with  a  full-wave  tomographic  imaging  effort,  as  good  earthquake  source 
mechanisms  are  essential  for  waveform-based  tomographic  inversions. 

6.2  An  Example  of  SGT-Based  Moment  Tensor  Solutions 

The  October  27,  2001  (Ms  5.5)  Yunnan  earthquake  was  recorded  by  a  temporary 
seismic  network  (de  la  Torre  et  ah,  2007)  and  the  permanent  stations  at  regional  and 


Approved  for  public  release;  distribution  is  unlimited. 

13 


teleseismic  distances  (Fig.  14g).  Body  and  surface  waves  at  both  regional  and 
teleseismie  distanees  in  all  frequeney  bands  provide  useful  eonstraints  on  the  souree 
parameters,  though  the  signals  are  best  in  the  30-60  s  period  range  (Fig.  14a-d). 
Approximately  l/3^‘*  of  the  frequency-time  windows  are  rejected  in  the  final  inversion 
because  of  their  low  signal-to-noise  ratios  (SNR)  or  correlation  coefficients.  The  souree 
loeations  in  the  EHB  and  GCMT  eatalogs  differ  by  '-'50  km  (Fig.  14f).  For  moderate  and 
small  earthquakes  like  this  one,  the  differenee  in  location  is  unlikely  the  effects  of  finite 
spatial  rupture,  but  dominated  by  location  errors  (Smith  &  Ekstrom,  1997).  Our  new 
solution  is  southeast  of  the  EHB  loeation  and  southwest  of  the  GCMT  solution,  though 
the  EHB  location  is  within  the  one-a  uncertainty  region  defined  by  waveform  misfit  and 
both  the  EHB  and  GCMT  loeations  are  within  the  95%  confidenee  limit.  The  depth  of 
the  preferred  solution  (8±3  km)  is  significantly  shallower  (at  the  1-g  uneertainty  level) 
than  the  GCMT  depth  (15  km)  and  the  EHB  depth  (16.7±3.2  km).  The  waveform  fit  at 
the  GCMT  and  EHB  depths  (15-17  km)  is  substantially  worse  than  the  fit  at  the  best 
solution  (Eig.  14e).  Compared  to  the  traveltime-based  EHB  solution  and  long-period 
waveform-based  GCMT  solution,  the  depth  resolution  in  this  study  comes  from 
variations  in  waveform  as  a  function  of  the  event  depth  and  shorter-period  body  waves 
used  in  this  study.  If  we  use  the  GCMT  moment  tensor,  the  waveform  misfit  is  nearly 
three  times  larger  than  that  of  the  best  SGT  solution  (1.64  of  GCMT  versus  0.52  of  SGT). 

Eigure  15a  shows  the  lowest  misfit  solution  at  eaeh  grid-seareh  depth.  The  solutions 
at  very  shallow  (2  km)  and  large  (>  18  km)  depths  are  noticeably  different  from  the 
global  best  solution.  In  contrast,  a  --lO  km  lateral  shift  in  loeation  does  not  affect  the 
moment  tensor  solution  signifieantly  (Eig.  15b).  This  low  sensitivity  to  lateral  shift  in 
loeation  is  partly  due  to  time  shift  in  waveform  fitting  to  account  for  errors  in  the  3D 
reference  model. 

A  close  comparison  of  the  observed  waveforms  and  the  synthetics  predicted  by  the 
GCMT  and  our  SGT  solutions  reveals  that  the  SGT  location  yields  waveforms  that  are 
better  matched  in  time.  To  illustrate  this,  we  zoom  in  on  the  first  arrivals  at  station  ESA 
and  CHTO,  and  the  teleseismic  S  wave  at  station  MDJ  (Eig.  16).  At  station  ESA,  the 
GCMT  and  EDSGT  syntheties  essentially  overlap  and  both  are  -'4±1  s  later  than  the 
observed.  The  delay  may  result  from  a  -'3%  lower- than-real  upper  mantle  P  veloeity  in 
the  referenee  model.  Since  a  significant  portion  of  the  event  to  ESA  path  cuts  through 
the  Indian  mantle  indentor,  the  veloeity  anomaly  seems  possible.  At  station  CHTO,  the 
onset  of  the  first  arrival  from  the  EDSGT  solution  matches  fairly  well  with  the  observed. 
Changing  the  source  to  the  GCMT  loeation  (and  the  GCMT  moment  tensor),  we  have  an 
onset  of  the  synthetie  first  arrival  that  appears  --lO  s  later  than  the  observed.  To  fit  the 
arrival,  the  upper  mantle  P  veloeity  along  the  event  to  CHTO  path  has  to  be  8%  higher 
than  in  the  reference  model,  a  value  too  large  to  be  likely  true.  At  station  MDJ,  the 
teleseismic  P  wave  is  below  the  SNR  threshold.  The  synthetic  teleseismic  S  wave  from 
the  SGT  solution  arrives  slightly  later  than  that  of  the  GCMT  solution  and  fits  the 
observed  better  in  time. 
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Figure  14.  Solution  of  the  October  27,  2001  (Ms  5.5)  Yunnan  earthquake,  (a-d)  Bandpass- 
filtered  observed  waveforms  (thick  gray  lines)  are  compared  with  the  synthetics  for  the  best 
source  solution  (red  lines).  The  frequency  bands  are,  from  left  to  right,  0.00833-0.0167,  0.0167- 
0.0333,  0.0333-0.0667,  and  0.0667-0.167  Hz.  Dark  thick  lines  are  time-shifted  observed 
waveforms  used  in  the  final  inversion.  For  the  first  arrivals  (in  c  and  d),  the  time  is  reduced  to 
align  the  arrivals  within  the  time  window.  The  amplitude  is  normalized  by  individual  trace,  (e) 
Normalized  waveform  misfit  as  a  function  of  depth.  Crosses  mark  the  grid-search  depths.  The 
shaded  region  indicates  the  1  -a  confidence  limit,  (f)  Contours  of  various  confidence  limits  at  the 
depth  of  the  global  best  solution  (8  km  for  this  event).  The  cross  marks  the  FDSGT  location,  the 
square  the  EHB  location,  and  the  diamond  the  GCMT  location.  From  the  innermost  contour 
surrounding  the  FDSGT  solution,  the  contours  are  the  68%  (1-g,  dark  line),  80%,  90%,  95%  (2- 
a,  dark  line),  and  99%  confidence  limits,  (g)  Triangles  and  cross  mark  the  stations  used  in  the 
source  inversion  and  the  event,  respectively.  Gray  contours  follow  the  sea  level,  and  2000  and 
4000  m  elevation. 
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Figure  15.  Variations  of  the  best  moment 
tensor  solution  as  a  function  of  depth  (upper 
panel)  and  lateral  locations  at  the  global  best 
solution  depth  (lower  panel).  The  global  best 
solution  is  shown  as  the  blue-and-white  beach 
ball.  The  GCMT  solution  (red-and-white 
beach  ball)  is  drawn  for  comparison.  The 
background  is  shaded  topographic  relief  with 
known  fault  traces  marked  by  dashed  lines. 


Figure  16.  A  comparison  of  the  observed 
waveforms  (thick  gray  lines)  and  the  synthetics  for 
the  SGT  (solid)  and  GCMT  (dashed)  solutions  at 
stations  LSA,  CFITO  and  MDJ  (triangles).  The  200 
random  realizations  of  the  FDSGT  solution  are 
plotted  below  the  observed  waveforms  to  show  the 
range  of  waveform  variation  due  to  uncertainties  in 
the  moment  tensor  solution.  The  waveforms  have 
been  filtered  between  0.01667-0.0333  FIz,  but  not 
shifted  in  time.  The  cross  and  circle  near  26° 
latitude  and  101°  longitude  are  the  SGT  and 
GCMT  locations,  respectively.  The  channel  and 
time  window  (after  the  earthquake  origin)  are 
marked  above  the  waveforms. 


7.  RESULTS  AND  CONCLUSIONS 

We  have  developed  a  novel  framework  in  seismic  tomography  to  construct  velocity 
models.  This  new  framework  results  from  the  perfect  marriage  of  the 
multigrid/multilevel  method  and  full-wave  tomography  and  provides  an  efficient  way  to 
cumulatively  and  progressively  improve  the  Earth  model  to  as  high  resolution  as 
observational  constraints  warrant.  We  apply  this  framework  to  build  a  3D  velocity  model 
for  the  eastern  hemisphere  that  covers  Eurasia,  Africa,  Mid  East,  Australia,  and  the 
Indian  Ocean  (latitude  55°S-55°N;  longitude  30°W-158°E). 

Using  a  frequency-time  data  normalization  procedure,  we  show  that  it  is  possible  to 
extract  up  to  300-s  period  Rayleigh  waves  from  ambient  noise  recorded  by  temporary 
broadband  seismic  stations  and  up  to  600-s  period  signals  from  permanent  stations  with 
very  broadband  sensors  (Shen  et  al,  2012).  These  very  broadband  signals  help  extend 
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surface  wave  constraints  to  all  depths  of  the  crust  and  upper  mantle.  We  eolleeted  and 
proeessed  up  to  23  years  (1990-2012)  of  eontinuous  seismie  data  from  broadband  seismie 
stations  in  the  eastern  hemisphere  (latitude  55°S-55°N;  longitude  30°W-158°E). 

We  have  developed  a  full-wave  ambient  noise  tomographie  method  (Shen  and  Zhang, 
2010;  Gao  and  Shen,  2014).  Three  levels  of  finite-differenee  eomputation  and  inversion 
have  been  earried  out  at  the  hemispherieal,  eontinental,  and  regional  scales  to 
progressively  eonstruct  a  hierarehieal,  multi-resolution  P  and  S  veloeity  model  for  the 
eastern  hemisphere.  The  results  show  much-improved  resolution  compared  to  those  of 
eonventional  ambient  noise  tomography  and  earthquake  surfaee  wave  tomography 
(Figures  9  and  13-1 1).  The  new  model  provides  eonstraints  on  several  geologieal 
proeesses  in  the  region,  sueh  as  the  origin  of  the  East  Afriean  Rift,  the  extent  of 
underthrustisng  of  the  Indian  lithosphere  beneath  Tibet,  and  reeycling  of  water  by 
subdueted  slabs.  Detailed  diseussion  of  the  geologieal  implieations  of  the  model  is 
beyond  the  seope  of  this  report. 

We  extend  the  SGT-based  souree  moment  inversion  from  loeal  to  regional  seales  and 
obtain  moment-tensor  determinations  of  earthquakes  between  2001  and  2004  in  the 
southeast  Tibetan  plateau  (Shen  et  ah,  2014).  The  SGT  inversion  takes  advantage  of  the 
souree-reeeiver  reeiproeity  and  is  well  suited  for  (near)  real  time  earthquake  souree 
studies.  The  SGT  database  approaeh  also  makes  it  eomputationally  effieient  to  determine 
the  best  solution  in  a  grid  seareh  surrounding  the  referenee  loeation,  by  minimizing  the 
waveform  misfit  of  time-windowed  first  arrivals,  teleseismie  S  waves  and  surfaee  waves 
reeorded  at  loeal  to  teleseismie  distanees.  The  new  solutions  are  in  general  agreement 
with  the  global  eentroid  moment  tensor  solutions,  though  the  two  souree  loeations  may 
differ  by  10s  of  km.  Inspeetion  of  waveform  fit  in  detail  suggests  that  regional  stations, 
the  3D  referenee  model  and  short-period  body  waves  tend  to  provide  tighter  eonstraints 
on  the  loeations  of  the  earthquakes.  The  lowest  magnitude  of  the  earthquakes  in  this 
study  is  Mw  4.35.  As  the  erustal  veloeity  model  improves,  short-period  (<15  s)  surfaee 
waves  beeome  useful  and  the  SGT  method  ean  be  applied  to  smaller  earthquakes  at  local 
and  regional  seales. 

8.  RECOMMENDATIONS 

The  multigrid  full-wave  tomography  framework  provides  a  powerful  means  to 
progressively  add  multiple  seismie  eonstraints  (body  and  surfaee  waves,  reeeiver 
funetions,  refieetion/refraetion  profiles)  at  various  plaees  at  any  time.  To  build  on  this 
framework,  we  recommend  the  following: 

(1)  Develop  efficient  wave  simulation  methods  and  eodes  that  take  advantage  of 
reeent  advanees  in  eomputing  teehnology  (e.g.,  GPEls  and  many  eores); 

(2)  Provide  a  shared-use  eomputing  faeility  at  an  aeademie  setting  (e.g.,  a  university) 
for  researchers  in  the  eommunity  to  earry  out  full-wave  tomography  and 
exehange/integrate  models; 

(3)  Provide  broad  training  on  the  use  of  the  eurrent  full-wave  tomographie  methods; 

(4)  Require  that  all  new  models  are  validated  against  independent  seismie 
observations. 

(5)  Colleet  and  maintain  a  waveform  database  that  eontains  independently  ealibrated 
ground-truth  seismie  events  for  model  validation 
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